!> date: 2022-10-23
!> 读取 pif-out.h5part 文件
module sph_load_h5part

    use h5part_m
    use seakeeping_kinds, only: rk
    use seakeeping_filesystem, only: is_exist
    use seakeeping_error_handling, only: file_not_found_error, file_parse_error
    use sph_region_type, only: region_t
    implicit none

contains

    !> 读取 pif-out.h5part 文件
    !> @note 
    !> 1. 步骤：打开并解析 nml 文件，打开 h5part 文件，读取步数，依次读取每一步的数据
    !> 2. 未读取加速度等额外信息
    subroutine load_h5part(file, nml, vtk, region, func)
        use seakeeping_math, only: arange
        use sph_io_env, only: idx, ones
        character(*), intent(in) :: file        !! 文件名
        character(*), intent(in) :: nml         !! Namelist 文件
        character(*), intent(in) :: vtk         !! VTK 所在文件夹
        type(region_t), intent(out) :: region   !! 计算域
        external :: func                        !! 输出数据的隐式函数
        integer(8) :: fid, stat, ik, nstep

        if (.not. is_exist(file)) call file_not_found_error(file)
        fid = h5pt_openr(file)
        if (fid < 0) call file_parse_error(file, 'file open failed')

        if (.not. is_exist(nml)) call file_not_found_error(nml)
        call load_pif_nml(nml, region%nreal, region%nvirt, region%ntotal, region%dim, region%hsml)

        allocate (idx(0:region%ntotal), source=arange(start=0, end=region%ntotal))
        allocate (ones(region%ntotal), source=1_1)

        allocate (region%loc(region%ntotal, 3), &
                  region%vel(region%ntotal, 3), source=0.0_rk)
        allocate (region%mass(region%ntotal), &
                  region%rho(region%ntotal), &
                  region%p(region%ntotal), &
                  region%itype(region%ntotal))

        nstep = h5pt_getnsteps(fid)
        do ik = 1_8, nstep, 1_8
            stat = h5pt_setstep(fid, ik)
            stat = h5pt_setnpoints(fid, int(region%ntotal, 8))
            stat = h5pt_readdata_r8(fid, 'x', region%loc(:, 1))
            stat = h5pt_readdata_r8(fid, 'y', region%loc(:, 2))

            stat = h5pt_readdata_r8(fid, 'Vx', region%vel(:, 1))
            stat = h5pt_readdata_r8(fid, 'Vy', region%vel(:, 2))

            stat = h5pt_readdata_r8(fid, 'Mass', region%mass)
            stat = h5pt_readdata_r8(fid, 'Density', region%rho)
            stat = h5pt_readdata_r8(fid, 'P', region%p)
            stat = h5pt_readdata_i4(fid, 'Itype', region%itype)
            call func(vtk, int(ik), region)
        end do

        deallocate (region%loc, region%vel, region%mass, region%rho, region%p, region%itype)
        stat = h5pt_close(fid)

    end subroutine load_h5part

    !> Load particle information from pif.nml file <br>
    !> 读取 pif.nml 粒子信息
    impure subroutine load_pif_nml(nml_file, nreal, nvirt, ntotal, dim, hsml)
        use seakeeping_string, only: to_string
        character(*), intent(in) :: nml_file                !! path of pif.nml file <br>
                                                            !! 文件路径
        integer, intent(out) :: nreal, nvirt, ntotal, dim   !! Integer variables <br>
                                                            !! 整数变量
        real(rk), intent(out) :: hsml                       !! Smoothing length <br>
                                                            !! 光滑长度
        integer :: iunit, precision
        character(6) :: creator
        character(23) :: create_date

        namelist /pif_attr/ creator, create_date, precision, nreal, nvirt, ntotal, dim, hsml

        open (newunit=iunit, file=nml_file, status="old")
        read (iunit, nml=pif_attr)
        close (iunit)
        if (dim /= 2) call file_parse_error(nml_file, &
            &"To-SPH-2D 仅处理二维问题，h5part 文件粒子维度为："//to_string(dim))

    end subroutine load_pif_nml

end module sph_load_h5part
